Make some waves.

Lets pretend we sample some really ideal sound at a rate of 10 measurments per second. Let x be your time line in seconds.

example.rate = 10
x <- seq(-30, 70, by = 1/example.rate)
wav.freq = c(0.8, 1.5, 2)
wav.1 = sin(wav.freq[1] * 2 * pi * x)
wav.2 = sin(wav.freq[2] * 2 * pi * x)
wav.3 = rep(0, length(x))
xmid = seq(0, 20, by = 1/example.rate)
someMidXs = which(x>0 & x < max(x)/3)
xmid = x[someMidXs]
wav.3[someMidXs] = 1.6 * sin(wav.freq[3] * 2 * pi * xmid)
wav.sum = wav.1 + wav.2 + wav.3
waves = list("wave1" = wav.1,
                         "wave2" = wav.2,
                         "wave3" = wav.3,
                         "sum" = wav.sum)

The Forier transfrom is beautifully explained by 3blue1brown in his video “But what is the Fourier Transform? A visual introduction.” I’m using that as my bases for making the following function.

The fourier transformation

fourtrans = function(winding, g, t){
    # winding - scaler, the "winding frequency"
    # g - a vector describing the original wave.
    # t - a vector of times corresponding to g; g and t must be the same length
    i = complex(imaginary=1)
    mass = g * exp( -2 * pi * i * winding * t)
    return(sum(sum(mass))) # the double sum is diliberate.
}

Do the fourier transform on the example waves.

ex.freqs = seq(0.001, 3, 0.001)
sum.ft = sapply(ex.freqs, fourtrans, g=wav.sum, t=x)
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
for (i in 1:length(wav.freq)){
    abline(v=wav.freq[i], col=wav.colors[i])
    text(x=wav.freq[i], y=par("usr")[4], xpd=T, pos=3, col=wav.colors[i], labels=wav.freq[i])
}

I see what I want. I see peaks that corresond to my original frequencies. But how to extract them when I don’t know the values and I don’t know how many there are? If I take every peak, I will have too many.

Pick out peaks

getMaxima = function(vector){
    n=length(vector)
    left = vector[1:(n-2)]
    main = vector[2:(n-1)]
    right = vector [3:n]
    indecies = which(main > left & main > right)
    # these indeces are for 'main' which is shifted from the original by 1.
    return(indecies + 1)
}
rawMaxima = getMaxima(abs(sum.ft))

Using the raw values and taking EVERY maximum gives me 227 frequencies.

plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col="red", pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
legend(legend=length(rawMaxima), x="topright", text.col="red", pch=16, col="red")

I could order them by their intensity (y-axis), but without knowing how many there are, I might stop with only 2 or take on several too many.

Optimize smoothing parameter - span

I know I can count on having a lot of little hills that don’t matter, so I’m going take the 75th percentile and use that a bench mark for “surely small”. Any maxima that are no higher than that, can be considered noise.

small = quantile(abs(sum.ft)[rawMaxima], .75)

dotColor = rep("red", length(rawMaxima))
dotColor[abs(sum.ft)[rawMaxima] <= small] = "lightblue"

plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col=dotColor, pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
abline(h=small, lty=2, col="lightblue")
legend(legend=sum(dotColor=="red"), x="topright", text.col="red", pch=16, col="red")

Smooth out the values to remove unimportant peaks.

n.panels = 8
par(mfrow=c(n.panels,1), las=1, mar=c(0,3,0,0), oma=c(4,0,3,0))

plotSmoothMaxes = function(x, y, dotsAt, thresh){
    #x - the x values of points to plot
    #y - the y values of points to plot
    #dotsAt - the indices of points to highlight (the maxes)
    #thresh - the threshold y-value below which to ignore the maxes
    plot(x, y, type="l", las=1, xlab="frequency", xaxt="n")
    polygon(x=c(max(x), min(x), x), c(0,0, y), col="gray")
    dotCol = rep("red", length(dotsAt))
    dotCol[y[dotsAt] <= thresh] = "lightblue"
    abline(h=thresh, lty=2, col="lightblue")
    points(col=dotCol, pch=16, x=x[dotsAt], y[dotsAt])
}

span = seq(0.01, .2, length.out = 40)
n.span = length(span)
drawThese = round(seq(1, n.span, length.out = n.panels))
drawList = list()
numMaxes = rep(NA, length(span))
for (i in 1:n.span){
    ft.lo = loess(abs(sum.ft) ~ ex.freqs, span=span[i])
    smoother = predict(ft.lo, ex.freqs)
    maxes = getMaxima(smoother)
    numMaxes[i] = sum(smoother[maxes] > small)
    if (i %in% drawThese) {
        plotSmoothMaxes(x=ex.freqs, y=smoother, dotsAt = maxes, thresh=small)
        abline(v=wav.freq, col=wav.colors)
        legend(x="topright", legend=paste("span =", span[i]), bty="n")
    }
}
title("smoothed data with different 'span'", outer = T)

Using a span of 0.07 is just enough to only get the main peaks. We still get the right number of peaks even with a span of more than twice that. As we increase the span, we see the values we would get from our maxima shift away from the accurate values (the red dots move away from the red, green and gold lines).

plot(span, numMaxes, type="l", las=1)
abline(h=3, lty=2, col="red")

Now we need a rule a computer can follow that will help it come to similar conclusions about optimizing the span for a given clip. There is probably some really good algorithm for doing that, but for now this will do:

# What is the most common number of peaks (across different span parameters)
# That must be right.  What is the lowest span that gives me that many? That must be best.
getBestSpan = function(numMaxes, span, showplot=F){
    tb = table(numMaxes)
    tb = tb[order(tb, decreasing = T)]
    mode = as.numeric(names(tb)[1])
    w = which(numMaxes == mode)
    lowSpan = span[min(w)]
    if (showplot){
        plot(span, numMaxes, type = "l", las=1)
        abline(v=lowSpan, col="red")
        abline(h=mode, col="blue", lty=2)
        legend(x="topright", legend=paste("span =", lowSpan), bty="n")
        title("Optimal Span")
    }
    return(lowSpan)
}
lowSpan = getBestSpan(numMaxes, span)

Based on this logic, 0.068 is the optimum span parameter for this clip.

getBestFreqs = function(ft, freqs, lowSpan, small, showplot=F){
    ft.lo.best = loess(abs(ft) ~ freqs, span=lowSpan)
    smoother = predict(ft.lo.best, freqs)
    maxes = getMaxima(smoother)
    maxes = maxes[smoother[maxes] > small]
    freq.of.interest = freqs[maxes]
    if(showplot){
        plotSmoothMaxes(x=freqs, y=smoother, dotsAt = maxes, thresh=small)
    }
    return(freq.of.interest)
}
freq.of.interest = getBestFreqs(sum.ft, ex.freqs, lowSpan, small=small)

And based on that value for span, the frequencies of interest are:

## [1] 0.798 1.500 1.997

sliding window fourier

The peak for the wave 3 frequency is shorter and wider than the other two. This makes sense since it was only present in a small portion of the clip while the other two were present consitently. Lets apply a sliding window approach.

Define the windows.

winSize = 100 # how may values wide each window is
shift = 10 # how many values from the start of one window to the start of the next
len = length(wav.sum) # how many total values in the clip
winStarts = seq(1, len-winSize, shift)
winEnds = seq(winSize, len, shift)
nw = length(winStarts) # number of windows

For each window, do the fourier transform and save the value for each frequency.

#ex.freqs = seq(0.001, 1, 0.001)
shiftingFT = matrix(NA, ncol = length(ex.freqs), nrow=nw,
                                        dimnames=list(startTime=x[winStarts], frequency=ex.freqs))
for (i in 1:nw){
    win = wav.sum[winStarts[i]:winEnds[i]]
    shiftingFT[i,] = sapply(ex.freqs, fourtrans, g=win, t=x[winStarts[i]:winEnds[i]])
}
shiftingFT.real = abs(shiftingFT)

frequency intensity over time

midTime=round((x[winStarts]+x[winEnds])/2)
image(x=midTime, y=ex.freqs, z=shiftingFT.real, las=1, col=heat.colors(length(shiftingFT.real)))

intensity over time line plot

For just the frequences of intersest, plot their intensity over time (across the sliding windows).

freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(freq.of.interest)
plot(1,1, xlim=range(x), ylim=c(0,300), type="n",
         xlab="time", ylab="frequency intensity", las=1)
for (i in 1:3){
    fq = freqName[i]
    lines(x=midTime, y=shiftingFT.real[,fq], col=wav.colors[i])
}
legend(x="topright", legend = round(freq.of.interest,3), col=wav.colors, lty=1)

The only thing we did to optimize the window size and shift size was… well nothing. I looked at the picture and adjusted it. I don’t know any good rules there, except hope that the same parameters work well over most cases.

Generate recovered waves

Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.

rebuilt = matrix(ncol=length(freq.of.interest),
                                 nrow=length(x),
                                 dimnames=list(x=x, freq=freqName))
for (i in 1:length(freq.of.interest)){
    fq = freqName[i]
    #.1 was the lowest span that I tried that didn't produce an error
    intensityPerWindow.lo = loess(shiftingFT.real[,fq] ~ midTime, span=0.1) 
    amplitudeFactor = predict(intensityPerWindow.lo, x)
    rebuilt[,i] = sin(freq.of.interest[i] * 2 * pi * x) * amplitudeFactor
}

Normalize the amplitude to roughly match the original.

maxWavSum = max(wav.sum)
intensitySumPerSec = rowSums(shiftingFT.real[,freqName])
maxIntesitySum = max(intensitySumPerSec)
ampFactor = maxWavSum / maxIntesitySum
rebuilt = rebuilt * ampFactor

Plot recovered waves

Not perfect… but not too bad. I think real sound data will have more cycles per unit, and I think that will mean an increase in resolution. It would also look cleaner to have used a cutoff rather than scaling the amplitude by the intensity, but I don’t think that is as true for real data as it is for this simulation.

Here’s the original for comparison:

getTimes = function(clipVals, rate=44100){
    #clipVals - a vector representing a sound clip.
    #rate - sample rate of the sound clip (probably 44,100)
    # value => a vector representing the time measurements for those values
    indecies = 1:length(clipVals)
    return(indecies/rate)
}
times = getTimes(arf, rate=sample.rate)

freq = seq(20, 20000, 40)

arf.ft = sapply(freq, fourtrans, g=arf, t=times)
arfMaxes = getMaxima(abs(arf.ft))

# plot(freq, abs(arf.ft), type="l", las=1, xlab="frequency")
# title("Fourier transformed dog bark")
# polygon(x=c(max(freq), min(freq), freq), c(0,0, abs(arf.ft)), col="gray")
# points(col="red", pch=16, x=freq[arfMaxes], y=abs(arf.ft)[arfMaxes])

plotSmoothMaxes(x=freq, y=abs(arf.ft), 
                                dotsAt = arfMaxes, 
                                thresh=quantile(abs(arf.ft)[getMaxima(abs(arf.ft))], .95))

# Earlier I was interested in a lot of plotting.  
# Here, just put the important bits into a function, no plots
getFreqOfInterest.old = function(ft, freqs,
                                                         small = quantile(abs(ft)[getMaxima(abs(ft))], .95), 
                                                         span = seq(0.01, .2, length.out = 40)){
    # ft = fourier transformed data
    # ex.freqs - frequences to iterate over
    # small = frequencies with intensities lower than this are discarded
    # span - a vector of spans to iterate through
    n.span = length(span)
    numMaxes = rep(NA, length(span))
    for (i in 1:n.span){
        tryCatch(expr = {
            numMaxes[i] = suppressWarnings(expr={
                length(getBestFreqs(ft, freqs, lowSpan=span[i], small=small, showplot = T))
                })
        }, finally = {
            #rm(ft.lo); rm(smoother); rm(maxes)
        })
        
    }
    lowSpan = getBestSpan(numMaxes, span, showplot = T)
    freq.of.interest = getBestFreqs(ft, freqs, lowSpan, small=small, showplot = T)
    return(freq.of.interest)
}
getPeakFreqs <- function(ft, small=NULL){
    if (is.null(small)){
        fracMax = max(ft)/10
        quant = quantile(abs(ft)[getMaxima(abs(ft))], .95)
        small = max(fracMax, quant)
    }
    return(findpeaks(abs(ft), minpeakheight = small, 
                                     nups = 2, ndowns = 2, minpeakdistance = 5))
}
freqPeaks = getPeakFreqs(ft=abs(arf.ft))
arfFOI = freq[freqPeaks[,2]]

#arfFOI = getFreqOfInterest(ft=arf.ft, freqs=freq)

arfFOI

## [1] 860 500

sliding window fourier

For each window, do the fourier transform and save the value for each frequency.

getShiftingFT = function(rawSound, sampleRate, 
                                                 winSize = 5000, shift = 1000,
                                                 freqs=freq){
    # rawSound - values from a sound file
    # sampleRate - number of values per second in rawSound
    # winSize - number of values per window
    # shift - number of values to shift over to start the next window
    # freq - freqencies to measure
    len = length(rawSound) # how many total values in the clip
    times = getTimes(rawSound, rate = sampleRate)
    winStarts = seq(1, len-winSize, shift)
    winEnds = seq(winSize, len, shift)
    midTime = (times[winStarts] + times[winEnds])/2
    nw = length(winStarts) # number of windows  
    
    shiftingFT = matrix(NA, ncol = length(freqs), nrow=nw,
                                            dimnames=list(midTime=round(midTime,3), frequency=freqs))
    
    for (i in 1:nw){
        win = rawSound[winStarts[i]:winEnds[i]]
        shiftingFT[i,] = sapply(freqs, fourtrans, g=win, t=times[winStarts[i]:winEnds[i]])
    }
    shiftingFT.real = abs(shiftingFT)
    return(shiftingFT.real)
}
arfMatrix = getShiftingFT(rawSound = arf, sampleRate = sample.rate)

frequency intensity over time

midTime = as.numeric(rownames(arfMatrix))
image(x=midTime, y=freq, z=arfMatrix, las=1, col=heat.colors(length(arfMatrix)))

Plotting all of the frequencies shows a lot of blank space, just show the freqencies up to 2x the highest freq of interest.

maxFreqShow = max(arfFOI) * 1.2
w = which(freq <= maxFreqShow)
plotMatrix = arfMatrix[,w]
image(x=midTime, y=as.numeric(colnames(plotMatrix)), z=plotMatrix, las=1, col=heat.colors(length(plotMatrix)), ylab="frequency (Hz)")

intensity over time line plot

#freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(arfFOI)
arf.col = wav.colors # they can just happen to be the same
plot(1,1, xlim=range(midTime), ylim=c(0,max(arfMatrix)), type="n",
         xlab="time", ylab="frequency intensity", las=1)
for (i in 1:length(arfFOI)){
    fn = freqName[i]
    lines(x=midTime, y=arfMatrix[,fn], col=arf.col[i])
}
legend(x="topright", legend = round(arfFOI,3), col=arf.col, lty=1)

Generate recovered waves

Hear each frequency alone.

for (i in 1:length(arfFOI)){
    toneTime = c(times, times+max(times), times+(2*max(times)))
    s = sin(arfFOI[i] * 2 * pi * toneTime)
    tone = audioSample(s, rate=sample.rate, bits=16)
    play(tone)
    wait(2)
}

Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.

rebuilt = matrix(ncol=length(arfFOI),
                                 nrow=length(times),
                                 dimnames=list(time=times, freq=freqName))
for (i in 1:length(arfFOI)){
    fq = freqName[i]
    intensityPerWindow.lo = loess(c(0,arfMatrix[,fq],0) ~ c(min(times),midTime,max(times))) 
    amplitudeFactor = predict(intensityPerWindow.lo, times)
    rebuilt[,i] = sin(arfFOI[i] * 2 * pi * times) * amplitudeFactor
}

Normalize the amplitude to roughly match the original.

maxWavSum = max(arf)
intensitySumPerSec = rowSums(arfMatrix[,freqName])
maxIntesitySum = max(intensitySumPerSec)
ampFactor = maxWavSum / maxIntesitySum
rebuilt = rebuilt * ampFactor

Combine generated frequencies into a sound

arfRebuilt = rowSums(rebuilt)
artArf = audioSample(arfRebuilt, rate=sample.rate, bits=16)

plot(start:end, arfRebuilt, type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")

play(artArf)

play(b1[start:end])
plot(start:end, b1[start:end], type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")

library("seewave")

## Warning: package 'seewave' was built under R version 3.4.4

tm = timer(b1,threshold=5,msmooth=c(50,0))

length(tm$s)

## [1] 264

This found 264 peaks, and I don’t see an easy way to merge adjascent peaks.

Smothing by maxes.

numB1 = abs(as.numeric(b1))
names(numB1) = getTimes(numB1)
nextTime = min(getTimes(numB1)) + max(getTimes(numB1))

vals = numB1
n = 9
naturalFloor = rep(NA, n)
inflatedShape = as.list(naturalFloor)

plot(x=getTimes(numB1), y=numB1, type="n", col="darkblue", las=1, xlim=c(1.1, 1.3))
for (i in 1:n){
    vals = vals[getMaxima(vals)]
    #tv = table(vals)
    #tv = tv[order(tv, decreasing = T)]
    naturalFloor[i] = min(vals)
    inflatedShape[[i]] = vals
    lines(x=as.numeric(names(vals)), y=vals, type="l", col=i)
    #hist(vals, main=i, 200)
}
legend(x="topright", lty=1, col=1:n, legend=1:n)

#plot(x=as.numeric(names(b1m1)), y=b1m1, type="l", col="darkblue", xlab="seconds", las=1)
nextTimeLab = as.character(nextTime)
inflatedShape = lapply(inflatedShape, function(x){x1=c(0, x, 0); names(x1)[1]="0"; names(x1)[length(x1)]=nextTimeLab; return(x1)})

Collect amplitude envelopes.

numB1 = abs(as.numeric(b1))
timesB1 = getTimes(numB1)

xlim=NULL

plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)

maxAmp = max(numB1)
minpeakheight = 0.1 * maxAmp
threshold = 0.01 * maxAmp
#threshold = naturalFloor[3]
allEnv = findpeaks(numB1, nups=5, ndowns=5, 
                                     minpeakheight = minpeakheight, 
                                     threshold = threshold, 
                                     minpeakdistance = sample.rate * .1)
allEnv = data.frame(allEnv)
names(allEnv) = c("height", "peak", "start", "end")
allEnv$peakTime = getTimes(numB1)[allEnv$peak]

points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")

xlim=c(1.10, 1.5)
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")

I wasn’t happy with the results from timer, and I didn’t like the output format. I was able to get reasonable peaks with findpeaks, but I didn’t like the start and end that was assigned to each peak. I want to use the peaks from findpeaks, but define the start and stop of each peak based on when the sound returns to some threshold.

To do this, I want to look for continuous areas of low values.

A value is only true if the surrounding n values are also true.

noLonelyTrue = function(bool, n=1, left=n, right=n, naEdges=NA){
    # bool - a vectore of boolean values
    # n - how many values on either side must be true for an original value to still be true
    # left, right - number of values to the left (preceeding) or right (after) that must be true
    end = length(bool)
    # matrices are faster than loops.  
    # any given row in this matrix can be read as "is the original value true, how about the one before it, and before that..."
    # the offset by one is acheived by adding a value at the end.
    toTheLeft = matrix(data=rep(c(bool,NA),left+1)[1:(end * (left+1))], nrow=length(bool), ncol=left+1, byrow = F)
    # same as above, but the other direction
    toTheRight = matrix(data=rep(c(rev(bool), NA), right+1)[1:(end * (right+1))], nrow=length(bool), ncol=right+1, byrow = F)
    
    # true values are treated as 1 in rowSums
    newBool = rowSums(toTheLeft)==left+1 & rev(rowSums(toTheRight))==right+1
    
    # the first and last n values don't have enough neighbors, by default they are NA values.
    if (!is.na(naEdges)){
        newBool[1:left] = naEdges
        newBool[(length(newBool)-right+1):length(newBool)] = naEdges
    }
    return(newBool)
}
bool = c(T,T,F,T,T,T,T,T,T,T,T,F,F,T,T,T,T)
noLonelyTrue(bool, n=1)

##  [1]    NA FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## [12] FALSE FALSE FALSE  TRUE  TRUE    NA

Now define start and stop of peask based on when there is a continuous low value.

plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")

# to test for values in pauses, use the level 3 smooth-by-maxima values
# shape.lo = loess(inflatedShape[[3]] ~ names(inflatedShape[[3]]))
# shape3 = predict(shape.lo, getTimes(numB1))
# threshold = max(threshold, min(shape3) * 5)
# lines(getTimes(numB1), shape3, col="orange", lwd=3, lty=2)

#which values are below threshold
peak=allEnv$peak
nrow=length(numB1)
ncol=length(peak)
bt1 = numB1 < threshold
bt2 = noLonelyTrue(bt1, n=20, naEdges = F)
bt = matrix(data=bt2, byrow = F, nrow=nrow, ncol=ncol)

beforePeak = rep(1:length(numB1), length(peak)) < rep(peak, each=length(numB1))
bp = matrix(data=beforePeak, byrow=F, nrow=nrow, ncol=ncol)
isBefore = bt & bp
startInd = apply(isBefore, 2, function(x){max(which(x))})
abline(v=getTimes(numB1)[startInd], col="green")
isAfter = bt & !bp
endInd = apply(isAfter, 2, function(x){min(which(x))})
abline(v=getTimes(numB1)[endInd], col="red")

xlim=NULL
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")
abline(v=getTimes(numB1)[startInd], col="green")
abline(v=getTimes(numB1)[endInd], col="red")